% Matlab code to estimate a flexible mixed logit model in wtp space.
% Written by Kenneth Train, first version Oct 21, 2015, 
% Heather Snell added discounting component April 29, 2019.
 
% Do not change the next line. It sets the global variables.

 format shortG
 clear all

% Sets the global variables.
global NP XMAT PROBS Z NZ NCS NDRAWS NV BETAS ZTYPE CrossCorr YesGPU COEF NROWS discountfun T modelname subset interaction

%Setting options and names.  Other options are not set up within this
%folder.
discountfun = "QuasiHyper";
subset = 0; %for full data set
estimate = 1;
interaction = 5; % 0 for no interactions, 
                 %1 for interaction with FullAtt, 
                 %2 for interaction with HighConseq
                 %3 for interaction with NotFullAtt
                 %4 for interaction with LowAtt
                 %5 for interaction with NoConseq
                 %6 for interaction with LowConseq
modelname = 'NoConseq';

diary off
delete NoConseq.out
diary NoConseq.out

filename = 'NoConseq.workspace.mat';
WantHessian=1;
WantBoot=0;

R_Range = [0.15 1;
           -0.5 0.5];
 P_Range = [0.2 1;
           -0.9 0.9]; %actually the beta range
 WTP_Range=[ 0  0.25; %price range
               -14	 12; %cap
                -18	 18; %sw
                -5	 6; %buffer
                -10.0187 10;%infra
                 0.6	 0.7174; %jobs
                -5.7674	 7.75; %quality
                -3.0512	 4.36;
                 -12	 10.25;
                -16	 16;
                -6	 6;
                -10.0187 10;
                -0.6	 0.7174;
                -5.7674	 7.75;
                -3.5512	 4.2668];


T=65; %max time periods;

CrossCorr=1;
NBins=20;

NReps=250;
NDRAWS=1000;
NGridPts=6001;

% ZTYPE=1 for polynomials.
%      =2 for step functions
%      =3 for splines
ZTYPE=3;

PolyOrder=8;
NLevels=6;  
NKnots=8;

% Load interaction data
XMAT = csvread('interaction_data.csv',1,0);

FirstRun = 1; %This is used and changed in Bootstapping, but don't change here
Nalts = 3; %number of alts per question
Nquest = 5; %number of questions each person answers
NROWS = size(XMAT, 1);
NP = NROWS/(Nalts*Nquest);
NCS = NP * Nquest; 

% 1. Person number (1-NP)            MUST BE THIS. DO NOT CHANGE.
% 2. Choice situation number (1-NCS) MUST BE THIS. DO NOT CHANGE.
% 3. Chosen alternative (1/0)        MUST BE THIS. DO NOT CHANGE.
% 4. Price 
% 5. CAP
% 6. SW
% 7. Buffer
% 8. Quality
% 9. Jobs
% 10. Infra
% 11. Wlife
% 12.-19 interaction with FullAttendance dummy
% 20-27 interactions with HighConseq dummy
% 28 FullAtt
% 29 HighConseq
% 30-37 interaction with LowAtt
% 38-45 interactions with NotFullAtt
% 46-53 interactions with NoConseq
% 54-61 interactions with LowConseq
% 62 LowAtt
% 63 NotFullAtt
% 64 NoConseq
% 65 LowConseq

if interaction==0
    if discountfun=="QuasiHyper"
        NAMES = {'Rate';  'Beta'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
                'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
    else
        NAMES = {'Rate'; 'PriceScale'; 'ASC.CAP'; 'ASC.SW'; 'Buffer';
        'Infra'; 'Jobs'; 'Quality'; 'Wlife'};
    end
elseif interaction==1
    intvec = csvread('fullatt.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.Att'; 'Beta'; 'Beta.Att'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.Att'; 'ASC.SW.Att'; 'Buffer.Att'; 'Infra.Att'; 'Jobs.Att'; 'Quality.Att'; 'Wlife.Att'};
        else
            NAMES = {'Rate';  'Rate.Att'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.Att'; 'ASC.SW.Att'; 'Buffer.Att'; 'Infra.Att'; 'Jobs.Att'; 'Quality.Att'; 'Wlife.Att'};
        end
elseif interaction==2
    intvec = csvread('highcon.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.Con'; 'Beta'; 'Beta.Con'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.Con'; 'ASC.SW.Con'; 'Buffer.Con'; 'Infra.Con'; 'Jobs.Con'; 'Quality.Con'; 'Wlife.Con'};
        else
            NAMES = {'Rate';  'Rate.Con'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.Con'; 'ASC.SW.Con'; 'Buffer.Con'; 'Infra.Con'; 'Jobs.Con'; 'Quality.Con'; 'Wlife.Con'};
        end
elseif interaction==3
    intvec = csvread('notfullatt.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.NotFullAtt'; 'Beta'; 'Beta.NotFullAtt'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.NotFullAtt'; 'ASC.SW.NotFullAtt'; 'Buffer.NotFullAtt'; 'Infra.NotFullAtt'; 'Jobs.NotFullAtt'; 'Quality.NotFullAtt'; 'Wlife.NotFullAtt'};
        else
            NAMES = {'Rate';  'Rate.NotFullAtt'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.NotFullAtt'; 'ASC.SW.NotFullAtt'; 'Buffer.NotFullAtt'; 'Infra.NotFullAtt'; 'Jobs.NotFullAtt'; 'Quality.NotFullAtt'; 'Wlife.NotFullAtt'};
        end
elseif interaction==4
    intvec = csvread('lowatt.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.LowAtt'; 'Beta'; 'Beta.LowAtt'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.LowAtt'; 'ASC.SW.LowAtt'; 'Buffer.LowAtt'; 'Infra.LowAtt'; 'Jobs.LowAtt'; 'Quality.LowAtt'; 'Wlife.LowAtt'};
        else
            NAMES = {'Rate';  'Rate.LowAtt'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                 'ASC.CAP.LowAtt'; 'ASC.SW.LowAtt'; 'Buffer.LowAtt'; 'Infra.LowAtt'; 'Jobs.LowAtt'; 'Quality.LowAtt'; 'Wlife.LowAtt'};
        end
elseif interaction==5
    intvec = csvread('nocon.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.NoCon'; 'Beta'; 'Beta.NoCon'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.NoCon'; 'ASC.SW.NoCon'; 'Buffer.NoCon'; 'Infra.NoCon'; 'Jobs.NoCon'; 'Quality.NoCon'; 'Wlife.NoCon'};
        else
            NAMES = {'Rate';  'Rate.NoCon'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.NoCon'; 'ASC.SW.NoCon'; 'Buffer.NoCon'; 'Infra.NoCon'; 'Jobs.NoCon'; 'Quality.NoCon'; 'Wlife.NoCon'};
        end
elseif interaction==6
    intvec = csvread('lowcon.csv'); %make sure no col names in file
        if discountfun=="QuasiHyper"
            NAMES = {'Rate';  'Rate.LowCon'; 'Beta'; 'Beta.LowCon'; 
                'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.LowCon'; 'ASC.SW.LowCon'; 'Buffer.LowCon'; 'Infra.LowCon'; 'Jobs.LowCon'; 'Quality.LowCon'; 'Wlife.LowCon'};
        else
            NAMES = {'Rate';  'Rate.LowCon'; 'PriceScale'; 
                'ASC.CAP'; 'ASC.SW'; 'Buffer'; 'Infra'; 'Jobs'; 'Quality'; 'Wlife';
                'ASC.CAP.LowCon'; 'ASC.SW.LowCon'; 'Buffer.LowCon'; 'Infra.LowCon'; 'Jobs.LowCon'; 'Quality.LowCon'; 'Wlife.LowCon'};
        end

end
        
IDPRICE=4;
% Get number of variables with flexible distributions
NV = size(NAMES, 1);
%Put supports together into COEF
COEF = [R_Range;P_Range;WTP_Range];

ThisSeed=1234;
rng(ThisSeed);

%Do not change the following lines. They calculate the number of Z variables
if ZTYPE==1 
    NZ=PolyOrder*NV;
    if CrossCorr==1
        if discountfun=="QuasiHyper" & interaction==0
            NZ=NZ+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
        elseif discountfun=="QuasiHyper" & interaction~=0
            NZ=NZ+(7)*(6)/2;
        else 
            NZ=NZ+(NV-2)*(NV-3)/2; %will not correlate price or rate
        end
    end
elseif ZTYPE==2
    NZ=(NLevels-1)*NV;
    if CrossCorr==1
            if discountfun=="QuasiHyper" & interaction==0
               NZ=NZ+2*(NV-3)+(NV-3)*(NV-4)/2; %will not correlate price or rate or beta
            elseif discountfun=="QuasiHyper" & interaction~=0
                NZ=NZ+2*(7)+(7)*(6)/2; %will not correlate price or rate or beta or interactions
            else
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
elseif ZTYPE==3
    NZ=(NKnots+1)*NV;
    if CrossCorr==1
            if discountfun=="QuasiHyper"
               if interaction~=0
                   NZ = NZ+2*(7)+(7)*(6)/2;
               else
                   NZ = NZ+2*(NV-3)+(NV-3)*(NV-4)/2;%will not correlate price or rate or beta
               end
            else 
               NZ=NZ+2*(NV-2)+(NV-2)*(NV-3)/2; %will not correlate price or rate
            end
    end
end

%Set the starting values for the coefficients of the Z variables
beta=zeros(NZ,1);  

YesGPU=1;

MAXITERS=5000;

%CONVERENCE TOLERANCE.
% Set =[] to use matlab defaults. default is 1e-6
PARAMTOL=[];
LLTOL=[];
tic;
if discountfun ~= "QuasiHyper"
    if interaction~=0 %interaction non-Qh
        CreateDrawsDisWtpProbsint;
        CheckRate;
        if check_ok==1
           CreateZint;
           disp(' ');
           disp(['Data setup took ' num2str(toc./60) ' minutes.']);
           disp(' ');
           if estimate==1
               DoEstimationRate;
               paramhat_orig = paramhat;
               Figure_OutputInt;
                if WantBoot==1
                    Boot
                end
            end
        end 
    else %interaction is 0 ---regular non-QH
        CreateDrawsDisWtpProbs;
        CheckRate;
        if check_ok==1
            CreateZ;
            disp(' ');
            disp(['Data setup took ' num2str(toc./60) ' minutes.']);
            disp(' ');
            if estimate==1
                DoEstimationRate;
                paramhat_orig = paramhat;
                Figure_Output9;
                if WantBoot==1
                    Boot
                end %end if boot
            end % end estimate
        end %end if check ok
    end % end interaction if else

else %discount is QH
    if interaction == 0
        CreateDrawsDisWtpProbsQH;
        CheckRate;
        if check_ok==1
            CreateZQH;
            disp(' ');
            disp(['Data setup took ' num2str(toc./60) ' minutes.']);
            disp(' ');
            if estimate==1
                DoEstimationRate;
                paramhat_orig = paramhat;
                Figure_Output10;
                if WantBoot==1
                    Boot
                end %end boot
            end
        end %end check ok

    else %interaction occurs
        CreateDrawsDisWtpProbsQHint;
        CheckRate;
        if check_ok==1
            CreateZQHint;
            disp(' ');
            disp(['Data setup took ' num2str(toc./60) ' minutes.']);
            disp(' ');
            if estimate ==1
                DoEstimationRate;
                paramhat_orig = paramhat;
                Figure_OutputInt;
                if WantBoot==1
                    Boot
                end %end boot
            end %end estimate
        end
    end
end
save(filename)
diary off
